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We investigate an open XX Z spin 1 /2 chain driven out of equilibrium by coupling with boundary reservoirs 
targeting different spin orientations in XY plane. Symmetries of the model are revealed which appear to be 
different for spin chains of odd and even sizes. As a result, spin current is found to alternate with chain length, 
ruling out the possibility of ballistic transport. Heat transport is switched off completely by virtue of another 
global symmetry. Further, we investigate the model numerically and analytically. At strong coupling, we 
find exact nonequilibrium steady state using a perturbation theory. The state is determined by solving secular 
conditions which guarantee self-consistency of the perturbative expansion. We find nontrivial dependence of 
the magnetization current on the spin chain anisotropy A in the critical region | A| < 1, and a phenomenon of 
tripling of the twisting angle along the chain for narrow lacunes of A. 



I. INTRODUCTION 



Manipulating of quantum systems which consist of a few quantum dots or quantum bits forms a basis of a functioning of any 
quantum computing device. Recent experimental advances allow to assemble and manipulate nanomagnets consisting of just 
a few atoms and perform measurements on nanostructures with an atomic resolution. However, a theoretical understanding of 
microscopic quantum systems out of equilibrium (e.g. under constant pumping or continuous measurement by a quantum probe) 
is far from being complete, apart from simplest cases like a single two-level system or a quantum harmonic oscillator under 
external pumping or in contact with a reservoir. 

Coupling of a quantum system with an environment (or with a measuring apparatus) is described under standard assumptions 
El, lU in the framework of a Lindblad Master equation (LME) for a reduced density matrix, where a unitary evolving part 
is complemented with a Lindblad dissipative action. Under the LME dynamics, a system with gradients evolves towards a 
nonequilibrium steady state (NESS) with currents. 

The Heisenberg model of interacting spins is an oldest many-body quantum model of a strongly coiTelated systems, important 
both from a theoretical and an experimental point of view. Energy and magnetization transport in the XX Z Heisenberg spin 
chain is extensively studied^by different methods, see B BI] for reviews, inclu ding Bethe Ansatz 10. bosonization M, 

J^, Lagrange multipliers tlQ], exact diagonalization lllLll2 1. Lanczos method ifisll . quantum Monte Carlo ifTill . DMRG ifisl 
[iql . Investigation of nonequilibrium spin chains from the point of view of Markovian Master equations contributes to the 
"equilibrium" studies. In quasi one-dimensional spin chain materials like SrCu02 many transport characteristics are measurable 
experimentally [fis!]. 

Investigation of NESS in driven spin chains is hindered by enormous technical difficulties due to an exponentially growing 
Hilbert space, although time-dependent DMRG method allows to reach somewhat larger sizes IitIi . In many studies the NESS 
are produced by coupling a spin chain at the edges to boundary baths of polarization aligning boundary spins along the anisotropy 
axis Z, so that the first and the last spin tend to be antiparallel [31.]. For this particular choice of boundary baths, a substantial 
progress in understanding the properties of a driven XX Z model have been recently achieved lfl9l - l24[ l26l - l30ll . including an 
exact solution of the problem in case of strong boundary driving ||23|, see also ll24ll and references therein. Note that within this 
boundary setup, one has the transport of a z-component of magnetization j^, caused by boundary gradients along the z-axis, 
thus probing diagonal elements of a diffusion matrix and a conductivity tensor. 

In present paper we carry out a systematic study of nonequilibrium steady states in driven XX Z spin chain resulting from 
boundary gradients in the XY direction, transverse to the anisotropy axis. In the context of spin and energy transport and for 
appropriate limits, this would correspond to probing the nondiagonal elements of the diffusion and conductivity tensors. To 
this end, we bring an XX Z spin chain in contact with reservoirs aligning the boundary spins along the X-axis at one end and 
along the F-axis in the other end, thus trying to impose a twisting angle in XF-plane of 7r/2 between the first and the last spin. 
Within the LME formalism, twisting angles a = ot a = tt result in a complete suppression of the magnetization current, see 
a remark after Eq.(fT6]l. making a setup with a = tt/2 the simplest nontrivial one. We derive general symmetries of the model 
which indicate that magnetization current should alternate its sign with system size, for any value of boundary coupling. On this 
basis we rule out the possi bility of a ballistic spin current, which is typically observed in integrable XX Z model under other 
kinds of perturbations jzsll . lSSll . We argue, further, that the phenomenon of the spin current alternation is rather general and in 
particularly is not related to the integrability. 

To confirm further our findings we investigate the nonequilibrium steady state of the XX Z model numerically and analyti- 
cally, for varying amplitudes of boundary couplings and boundary gradients. We find complete agreement with our predictions 
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Figure 1: Schematic layout of the problem. The chain of two-level systems coupled at the boundaries to the reservoirs. Dissipation introduced 
by the reservoirs is described by a quantum Master equation (TJ 

which in particular are confirmed analytically, in the regime of strong coupling. In this regime, full analytic solution of the 
Lindblad equation for large Lindblad amplitudes demonstrates a very nonmonotonic dependence of the magnetization current 
on the system size and anisotropy. 

The plan of the paper is the following. In sec|II]we introduce the model. Symmetries of the Lindblad master equation are 
obtained in sec|IIIl Using the symmetries general NESS properties for arbitrary system size and coupling are revealed, including 
an admissibility of spin and energy currents. Parity-depending LME symmetries lead to a phenomenon of magnetization current 
sign alternation with the system size. In sec. |IV]our predictions are tested further by calculating the NESS analytically and 
numerically for small system sizes. The cases of XX model and of XX Z model are treated separately in sec lTV Al and sec II V B| 
Limits of weak coupling and of weak driving force are briefly discussed at the end of sec lIV CI Appendix contain necessary 
technical details. In conclusion we summarize our findings. 

II. THE MODEL 

We study an open chain of N quantum spins in contact with boundary reservoirs. The time evolution of the reduced density 
matrix p is described by a quantum Master equation in the Lindblad form lit], 101 > ll29ll (here and below we set /i = 1) 



^^=-i[H,p]+T{Cl[p\+Cb[p\), (1) 

where H is the Hamiltonian of an extended quantum system, and Cl[p\ and Cft\p\ are Lindblad dissipators acting on spatial 
edges of the system. This setting, shown schematically in Figd] is a common starting point in studies of nonequilibrium 
transport, a particular choice of Cl and C r depending on the application. Often, a choice of the Lindblad action is operational, 
i.e. is determined by a requirement to favour a relaxation of a system or its part towards a target state with given properties, e.g. 
to a state with given polarization or temperature UtI] . In this way, one describes an effective coupling of the system, or a part 
of it, with a respective bath-environment. Note that within the quantum protocol of repeated interactions, ll29ll . lf34ll the equation 
([T]i appears to describe an exact time evolution provided that the coupling to reservoir is rescaled appropriately with the time 
interval between consecutive interactions of the system with the reservoirs. 

We specify the Hamiltonian to describe an open XX Z spin chain with an anisotropy A 

Af-l 

^ = E ("'k^Ui + <^M+i + AcT|a|+i) , (2) 
fc=i 

while Cl and Cfi are Lindblad dissipators favouring a relaxation of boundary spins k = 1 and k ~ N towards states described 
by one-site density matrices and pa, i.e. Cl[pl] = and Culpfi] = 0. General form of the matrices p^ and pfj is 
Pl,r = 5 + ctL.flC^ + Pl,r'^^ + 1l,rO'^^ where tr" are Pauli matrices, and a, 13, 7 are real constants satisfying Trp| ^ < 1 or 

^L,B. + f^L.R + Il.r — 1- ^ general choice of the Lindblad dissipator £^^7^ [.] = X^i j=i af^-'^ fo-'. cr^ — i{ . , cr-'cr*}), each 
of the 2x2 matrices p^ and pR, if not constant, belongs to at most one-dimensional manifold IBal . 

We choose boundary reservoirs which tend to align the spin at the left edge along the x direction and the spin at the right edge 
along the y direction. Consequently, the reservoirs try to establish a twisting angle of 7r/2 between the first and the last spin in 
the chain. Such a setting is achieved by taking the Lindblad action in the form 




r.(3M3H5H5M=H5>- 

1 2 N 
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W = - 2 ^ ^IW^] + E ^™/^^™, (3) 

m— 1 7n— 1 

1 ^ ^ 

m— 1 m— 1 

Here and are polarization targeting Lindblad operators, which act on the first and on the last spin respectively, 



W2 = W^K 



Vl^\^-T^{^l+^a^,). (4) 



In absence of the unitary term in ([T]) the boundary spins relax with a characteristic time T^^ to specific states described via 
one-site density matrices pL and pa, satisfying Cl[pl\ — and Cii[ppi\ = 0, where 

PL = \{I- ^^l) (5) 

PR^\{I + i^'<^n) ■ (6) 

From the definition of a one-site density matrix Pone-site = 5 + '^^at the Lindblad superoperators Cl and 

£ii indeed try to impose a twisting angle of tt/2 in XY plane between the first and the last spin Issll . The twisting gradient 
drives the system in a steady-state with currents. In the following we restrict to a symmetric choice 

< K = k' < 1. (7) 

The parameter k determines the amplitude of the gradient between the left and right boundary, and therefore plays a fundamen- 
tally different role from the coupling strength F. The limits k = 1 and k <C 1 will be referred to as the strong driving and weak 
driving case respectively. The two limits describe very different physical situations as exemplified in other NESS studies. 

The reason for choosing the Lindblad dissipators which impose a twisting gradient in XF-plane perpendicular to the Heisen- 
berg chain anisotropy is two-fold. Firstly, our choice is "orthogonal" and therefore complementary to usually considered bound- 
ary setup fl9ll - ll24ll . 112611 - lEsll . where alignment of spins at the boundary is parallel to the anisotropy axis Z. Our system is not 
mappable to a quadratic fermionic system, and is therefore not amenable to a analysis via a canonical quantization in the Fock 
space of operators lfl9ll . ll20ll . which can be used for usual boundary setup. Secondly, one realizes that imposing both X- and Y- 
boundary gradients is necessary to create a nonvanishing stationary Z-magnetization current. A boundary setup with parallel or 
antiparallel boundary spins, aligned in the same direction in XY- plane leads to zero steady magnetization current, see Remark 
after iT% . Consequently, by imposing a twisting angle in XF-plane of tt/2 we choose a minimal nontrivial setup. 

Transport properties in the XX Z spin chain H = hn,n+i are governed by spin and energy current operators which are 
defined by lattice continuity equations J^crf^ = jn^i^n - in,«+i, J^/i„,„+i = - Jn+i where 

i„,„. = 2(«„ - «), (8) 



We show next that the magnetization current j alternates its sign with the system size, while the energy transport is switched 
off completely, if k = k', despite the boundary gradients. Being locally conserved quantities, in the stationary state the currents 
(in,n+i) = ii {Jn) ~ ^'^ cqual across all the bonds. 

For the following is essential to note, that for our choice of the Lindblad dissipator ([3]), (|4|i, the system ([T]i has a unique 
nonequilibrium steady state. Its existence and uniqueness for any coupling F is guaranteed by the completeness of the algebra, 
generated by the set of operators { H, Vm , Wm , Kn 7 } under multiplication and addition ll39ll , and is verified straighforwardly 
as in H. 
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in. SYMMETRIES OF THE LINDBLAD EQUATION AND RESTRICTIONS ON ENERGY AND MAGNETIZATION 

CURRENTS 

Let us denote by p{N,A,t) a time-dependent solution of the Lindblad equation ([T]i, ^ for a system of N sites, and by 
p{N, A) ~ limf^oo p{N, A, t) its steady state solution, where A is the anisotropy of the XX Z Hamiltonian (|2]i. There exist 
transformations which map one LME solution to another Some of the transformations depend on the parity of A^. For even N 
, introducing U = ~ Yi ® '^n' we find that p{N, — A, t) = Up*{N, A, t)U is a solution of the same Lindblad equation 

n odd 

with the anisotropy —A. We denote by p* a complex conjugated matrix in the basis where <t^ is diagonal. However, if a 
nonequiUbrium steady state is unique, then a relation follows which relates the NESS for positive and negative A, 

p(7V,-A) = C/p*(iV, A)C/foriVeven. (10) 

For odd N we obtain another relation for NESS, introducing a unitary transformation = (cr^)®" 

p{N,~A) :^J:yUp*{N,A)UY.yfoiNodd (11) 

Another LME symmetry does not depend on parity of N and maps a NESS onto itself provided that an additional condition 
O is satisfied: 

p{N, A) = E.UrotRpiN, A)i?C/+tE,, (12) 

where R{A (g) B ... (g) C) = (C (X) .... A)R is a left-right reflection, and the diagonal matrix Urot = diag{l, i)®" is a 
rotation in XY plane, Urot^ynU^^i = cr^, UrotiJnU^ot ^ ^"'n' ™d = (cr^)®" . Note also that at the point A = symmetries 
( [Tol l and (fTTl i map a NESS onto itself, and are independent on (fT2l) . 

The relations ( fTOl ) - ( fTSI ) impose restrictions on the NESS density matrices and respectively on all observables including the 
magnetic and energy currents j = Tr{p]n,n+i) and J'^ = Tr{pJ^) where jn.n+i and are given by (H) and (|9]l. We drop 
the subscripts since in the steady state the currents j, are equal across all bonds. By virtue of the symmetry ( fTOl i. we obtain 
for N even 

j(-A)=TrOV„+ip(iV,-A)) = 

^Tr{3r,^n+iUp*{N,A)U) = 
= -Tr{U]^^n+iUp{N,A))* = 
= rrOV„+ip(iV,A)) =j(A). 

In the last passage we used (i) j* = — jn,n+i in the basis where <t^ is diagonal (ii) Ujn,n+iU = —jn.n+i (iii) the fact that 
an observable j is a real number. Analogously, for odd N, we obtain using (fTTT i. that j(— A) = — j(— A). So, we find that the 
magnetization current changes its parity from even to odd as a function of the anisotropy A: 

.7(A) = j(-A) for N even, (13) 
.7(A) = -j(-A) for iV odd. (14) 

We argue below, see sec lIV Al that at the free fermion point A = the spin current vanishes always, j{A = 0) = 0. For 
odd N it follows from ( fT4l i. while for even N it follows from a perturbative expansion. In addition, we make another suprising 
observation: due to (flJl i. (fl4l i the spin current alternates its sign with the system size since by increasing a system size by one 
unit — !■ iV + 1 the current dependence on A changes from even to odd or vice versa. For our case it amounts to 

sign j {A, N) = {-1)'^ for A < -1, (15) 

because in the region |A| < 1 the sign of j(A) can change depending on other parameters, see sec lIV Bl The alternating current 
effect is a consequence of an existence of different LME symmetries for even and odd number of sites ( fTOl i. (fTTI) and as such is 
quite general. Indeed the derivation of (fT3]),(fT4b does not depend neither on coupling F now on gradients k and k' as given by 
(|4]i. Also, the phenomenon of alternating current does not rely on an integrability of the XX Z model. Indeed, an addition of 
a staggered anisotropy to the Hamiltonian (|2]i breaks integrability but the Eqs. ( fTOl i. ( fTTI ) remain valid. Moreover, even though 
our conclusion relied on the validity of Eqs. ( fTSI ). (fT4l) . the alternating current was also observed in a recent study of a driven 
XX Z chain with an XY boundary setup different from (01), where one of two symmetries (EqdTSI)) was absent ll32ll . We suggest 
that the magnetization current sign alternation resulting from applying transverse gradients can be a rather general phenomenon, 
which might be observable under appropriate experimental conditions in artificially assembled nanomagnets consisting of just a 
few atoms Ii41i] . 
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Another consequence of the current sign alteration is that the spin current cannot be ballistic. Indeed, for large N we should not 
be able to differentiate between even and odd N, so that the stationary magnetization current j must vanish in the thermodynamic 
limit: 

jU^oo =OforallA,r. (16) 

This conclusion is a rather unexpected one, since in the XX Z integrable system, in the critical region | A| < 1 one usually 
expects a ballistic current llzsll . The magnetization current is also ballistic | A| < 1 for reservoirs creating boundary gradients 
along the anisotropy axis 1330 . 

Up to now we did not make use of the symmetry ( fT2] i. In fact, it does not impose any further restrictions for spin current, 
while the energy current changes its sign under the symmetry ( fT2] i and must therefore vanish 

J^^ ^Oforall Af,A,r. (17) 

Consequently, the energy and spin currents decouple: the spin current can flow, while the energy current is totally suppressed 
([TtI i. Vanishing of the energy current is a consequence of a symmetric choice k = k' dTji, and will be lifted if k 7^ n' . 

Remark. Note that if both boundary reservoirs were aligning boundary spins along the same direction in Xy-plane, the 
magnetization current will be completely suppressed for all A, despite the existence of boundary gradients. E.g. for a choice of 
the set of Lindblad operators Vi = a(<7^ + V2 = 0, Wi = (3{(7^ + icrf^), W2 = in ([3]i, which favours the alignment 
of the boundary spins along the X axis, the current suppression follows from a symmetry p — J^^P^x of the NESS, since 
Tr{jn^n+ip) = ~Tr{jn.n+i^xP^x)- Thus, our choice (O is the minimal nontrivial one, to observe a nonzero magnetization 
current. 

In the following we check our findings and analyze the model further by computing exact NESS numerically for finite F and 
analytically for F ^> co using a perturbative expansion. 



IV. STRONGLY DRIVEN XX Z CHAIN WITH XY TWIST: NUMERICAL AND ANALYTICAL STUDY. 

We performed studies of the LME for small system sizes N in order to support our findings (I13m7b . which should be valid 
for all N, A and F. For simplicity we restrict here to the strongly driven case k = k' = 1. Being a problem of exponential 
growing complexity, the steady state for a finite F is a very complicated function of the parameters F, A, for any iV > 3, and 
cannot reported otherwise than in a graphical form. We solve the linear system of equations which determine the full steady 
state, numerically, making use of the global symmetry ( fT2] i which decrease the number of uknowns by roughly the factor of 1/2. 

For iV = 3, 4 the magnetization currents as function of the anisotropy, for various values of F are reported in Fig l2l3l As 
expected from (fT3]l,(fT4li. the current is an odd (even) function of A for N ~ 3{N ~ 4). Interestingly, as the coupling F increases, 
the current j(F, A) approaches quickly a limiting curve drawn by a dashed line. For the dashed line, which corresponds to the 
the limit of infinitely strong couplings F — > cxd, a simple closed analytic form can be obtained. Here it is worthwhile to note that 
Lindblad Master equation dynamics ^ with large effective coupling to a reservoir F ;» 1 is realizable within the experimental 
protocol involving repeated periodic interactions of a boundary spin with a infinite bath system of identical ancilla ll29ll . 

j 



1.5 r 




-1.5^ 



Figure 2: NESS spin current j vs A for = 3 for different F — 0.5, 5, 10 (from bottom up). For finite F, the data are obtained numerically. 
Dashed line marks j{A) in the limit F — >■ cxj from i29i 
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Figure 3: NESS spin current j vs A for A'^ = 4 for different F — 0.5, 5, 50 (from bottom up). Tiie j{A) curves continue symmetrically in the 
negative A region. For finite F the data are obtained numerically. Dashed line marks j (A) in the limit F — >■ oo from OOt . 



To investigate the strong coupling limit F 1 we search for a stationary solution of the Lindblad equation in the form of a 
perturbative expansion 

PArsSs(A,r)=^ — j pfe(A). (18) 

As stressed in the end of Sec|lll the stationary solution pNEssi'^, T) is unique. The zeroth order term of the expansion (fTSb 
po = limr^oo linit^oo p(r, A, t) satisfies ClrIp'^*^] = (here and below we denote by Clr = Cl + Cji the sum of the 
Lindblad actions in This enforces a factorized form 



Po = PL ® t: + Mo{A) ® PR, (19) 




where pL and pn are one-site density matrices (|5]l, (|6]l and Mq is a matrix to be determined self-consistently later. We shall drop 
A- and A^-dependence in p for brevity of notations. We separate the identity matrix (|) for future convenience, so that 
Tr(po) = 1 and Mq is traceless. 

Inserting the ( fTSl l into dTJ, and comparing the orders of F"*^', we obtain recurrence relations 

i[H, pk] = ^jCLRPk+i , fc = 0, 1, 2... 

A formal solution of the above is pk+i = —"^^LniQk+i) where Qk+i ~ —i[H, pk]- Note however that the operator Clr has a 
nonempty kernel subspace, and is not invertible on the elements from it. The kernel subspace kcr(£Li?) consists of all matrices 
of type pi^ ® A ® PR where A is an arbitrary 2^^^ x 2^^^ matrix. Therefore a necessary and sufficient condition for Pk+i to 
exist is to require a null overlap 

[H, pfc] n ker^CLR) = (20) 

We shall name ( |20l i the secular conditions. For our choice of the Lindblad operator the secular condition are equivalent (see the 
Appendix) to the requirement of a null partial trace 

Tr^.,N{[H,Pk]) = Q, k = Q,l,2.... (21) 

The remaining missing ingredient of the perturbation theory is obtained by noting that pk+i is defined up to an arbitrary element 
from ker(£ijf), so we have 

Pk+i=2Cl\{i[H,pk\) + PL®Mk+i®PR, fc = 0,l,2... (22) 

Eqs ( fT9b . (I22I) and (I2TI) define a perturbation theory for the Lindblad equation ([T]). At each order of the expansion, we must 
satisfy the secular conditions ( |20| |. before proceeding to the next order. Construction of the inverse C^]^ is illustrated in the 
Appendix. 

Applying the perturbation theory to our model, we find that the unknown apriori matrices M2fe(A), Af2fe+i(A) are fully 
determined by secular conditions d^TT i for 2fc, 2k + 1, as will be illustrated further on an example. For any A 7^ the set of 
matrices {M;;(A)} is nontrivial. The case of zero anisotropy A = turns out to be special in that all = and is discussed 
separately below. 
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A. "Free fermion" point A = 0. Arbitrary A'^ 



For A = we can guess a zero approximation po for arbitrary N: 

2 



Po = PL ( ^ I ® PR- (23) 



Comparing with ( [19] ). we see that Mq — 0. All further — and all secular conditions ( |2TI ) are trivially satisfied at all orders. 
This fact was checked explicitly for small system sizes N < 7 and is conjectured to be true for arbitrary N. Consequently, the 
general recurrence ( |22] | is reduced to pk+i = 2££j^(j[i^, Pk]) for all k. At the first order of expansion, we obtain i[H, po] = Qi, 
where 



N-3 



/l = -cr (>$ (T (K) I - ) ()<)PR 



which is easily inverted pi = 2££]j,(Qi) = — Qi and so on. At all orders we notice Tr{jpk) = 0, which allows to conjecture 
that magnetization current vanishes for A = in spite of twisting boundary gradients, 

J I A=o = for all iV > 3 and all T. (24) 

Note that for odd N, the statement ( |24] | also follows from ( fT4l i. 

For iV = 3, the perturbative expansion ( fTSl ) can be summed up for all orders due to the fact that it closes at the third step: 
P3 = — 4pi. Consequently, ( fTSl l is rewritten as 

/I 1 \ °° 1 

P^=3,A=0(r) = PO + i^-pi + ^P2 j E (25) 

where Pi = -^{a^ ® a'' (g) pR + pL ® (gxj^) and p2 = icr^ ® / (g) pi? - icr^ » ct^ (g) cr"^ + icr"^ (g) cr"^ ® ti^ + jPl / » ct^. 
Summing up the series, we obtain 



PA^=3,A=o(r) = PO + + ^^^T^y (26) 



It can be verified straightforwardly that ( |26] | is a stationary solution of LME ([U for = 3, A = 0. The sum dZST l converges for 
r > 1. However, since the LME solution is unique, the expression (|26] | is valid for any F. In particular, for F — !• we obtain the 
weak coupUng Umit p^"^°-^[/S.) = limr-s-o linif^oo p(r, A, t), which is also unique by a continuity argument, 

P7V=3,A=o(r ^ 0) = p"'-^(A)|^^„ = PO + \p2. 

As expected, p""^"'^(A)|^u commutes with the XXO Hamitonian, [po + \p2, Hxxo] = 0, and corresponds to zero currents, 
j(F ^ 0) = 0. 

From the above example we see that perturbative series (fTsT l are useful even outside of their radius of convergence: if the 
perturbative series for pness or for a particular observable (/) ~ Tr[f p^ESs) can be summed, the result of the summation is 
valid for all F within its analyticity range, see also ( [35] ), ( [36l ). 



B. Exact NESS for small system sizes and arbitrary A. 

For any nonzero A, the matrices Mk from (l22T i are nontrivial and are determined from secular conditions (ISTT i. The secular 
conditions for two consecutive orders 2m, 2m + 1 reduce to a system of linear equations which determine and M2m+i- 
In such a way, the perturbative series ( fTSl l is uniquely defined at all orders. 

Below we present analytic results for A/o and small N which define the NESS through Eq.(fT9l). A/o and Mi are obtained 
by satisfying secular conditions ( 1211 1 for fc = and fc = 1. The procedure is illustrated below for the simplest case N = 3,4. 
For iV = 3, the most general form of the matrices Afo(A), Mi (A), by virtue of the symmetry ( fT2] i. is A/o = fio(f^ ~ o"^)- 
Ml = 6i((T^ — tr^) where 6o, &i are unknown constants. The secular conditions ( |2T] ) for fc = do not give any constraints on 
ho, while those for fc = 1 give two nontrivial relations 6i = and —A + (1 + 2A^)6o — from which both Mq and A/i are 
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determined. For TV = 4 each of the matrices Mq (A) , Mi (A) compatible with the symmetry ( fT2] l. contains 9 unknowns, all fixed 
by the secular conditions ( ITTT i for fc = 0, 1, and so on. For = 3, 4 we obtain 

Mo(A)|^^3^^^^(a--a^) (27) 

s(A)Mo(A)|^^4 - (^-^ + 2A4^ ® a- + Q _ 2A2^ 

- 4A2cr2^ ® cr^ + Q + 2A2^ {a'' ®I -I® ay) (28) 

A(4A2 + 5),^ , 
H ^ — -{I ®CT-^ ~(jy ®I) 

+ 2A((T^®CT^+cr2'®cry). 

where s(A) = 7 + 6A^ + 8A'*. Matrices i\/o obey the symmetries (fT0b-(fT2li. Steady state spin currents jo(A) = Tr{poj) in 
the hmit F oo are 



4A 

io(A)U^3 = Y^^, (29) 



•rA^I 4(1 + 4A2) 



7 + 6A2 + 8A4' 



while the energy current vanishes J'^ = as expected from ( [TtI i. We note that the current for iV = 4 is nonzero at A = 
which seems to contradict our prediction (|24] |. The contradiction arises because the limits F^^ — and A — s- do not commute. 
Indeed, we find that even orders of the perturbative expansion p2m contain terms which become singular for A — > 0. The leading 
singular term at the order 2m of the expansion is Tr{e'^"'' p2mj) ~ (e/A)^™, where e — F~^. Since the singularities become 
stronger with m, for any given infinitesimal e there will be a region | A| < e where the expansion dTSI l will diverge. The origin 
behind the non-commuting limits is the existence of an extra symmetry ( fTSl ) of the NESS at the point A = 0. Similar situation is 
encountered e.g. in Kubo linear response theory, where a nonergodicity at zero point in the momentum space if = is inflicted 
by an extra conservation law at this point ll42ll . Reversing order of the limits e — > and A — !> does yield j{0) = for any F, 
so that the current expression for = 4, F — !• oo has to be modified as 

Jo[^)\N=i~<^ 0, ifA = 0, 

while Mo|jY=4 A=o = as given bv (l23]l. For any finite F the dependence j( A) is smooth and does not contain any singularities, 
see Fig|3] We find the noncommutativity of the limits e — ^ 0, A — J> for all even system sizes N — 2m > 4, see also (|32] |. Since 
for large N one should not be able to differentiate between even and odd system sizes, we expect that the jo (A) for A ^ 1 
vanishes in the thermodynamic limit: limm^oo jo{^ ^ ^)\N=2m ^ 

The number of equations to solve ( l2Tl i at each order of perturbation theory grows exponentially with the system size N, 
although the symmetry ( fT2] i restricts the number of independent variables. The matrices Mq for > 5 are given by lengthy 
expressions and are not reproduced here. Instead, we report NESS currents at F ^> cxd, given for iV = 5, 6 by 



4A3(-3 + 12A2 + 23 A4 + 2^A^) 
1 - 8A2 - 6A4 + 120A6 + 96A8 - 26A4 + 2''A^^ 



JoIat^S ^;A4iir,nAfiin^;AS r>fiA4ir,7A12 ^^^^ 



„• I _ J 149-604A2+316A''-912A8+4992AS-5760Aio + 15360Ai2-2i'iAi*+2i3Ai6 ' ^ ^ 7^ U. 

Jo\n=6-[ 0, ifA-0 ^ 

The NESS currents jo for odd/even number of sites N are even/odd functions of A respectively, as expected from (l3[ . (fT4l) 
and depend rather nontrivially on the anisotropy A, see Fig|4] In particular for iV > 5 one observes sharp peaks in the 
jo (A) resembling resonances at which the current changes its sign and amplitude. One gets an insight by comparing av- 
erage magnetization profiles (ct;^) and (a^) for values of A inside and outside the resonance region. The average magne- 
tization at boundaries n ~ 1,A^ is fixed by the Lindblad reservoirs (erf) — !> 0, (crf^) 1 and {af) 1, (cr^) 
so that the lowest interpolating Fourier mode has the wavelength A/a — A{N — 1) where a is the lattice constant. In the 
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(b) 



Figure 4: NESS spin current vs A for = 5, 6 from Eqsl l31t ,ll32l) given by thick and thin line respectively (Panel (a)), and for A'' = 7 (Panel 
(b)), all for r 00. 

resonance region we observe that A' = A/3 becomes the dominant Fourier mode, see Fig|5jb). One can interpret this fre- 
quency tripUng as a triple increase of a twisting angle along the chain. Note, see Fig|5ta,b) that there is a 7r/2 phase shift 
between the dominant harmonics describing (cr^) and {a"!^) profiles. In particular, in the resonance region, Fig|3b), (<t^) 
K, — siii0(n — 1), (fT-() w — cos(/)(7i — 1), where = Stt/IO. In the meanfield approximation the respective magnetization 
current is jj^ip = '^{{<^n){<^n+i) ~ {'^n){'^n+i)) — — 2sin(/) k, —1.61803, which is amazingly close to the exact current value 
at the resonance jexact = —1.61857... With the increase of system size, one might expect an appearance of further resonances 
at wavelengths A' = A/5, A/7, ... corresponding to larger twisting angles ^ etc.. Why higher Fourier harmonics become 
stable under a strong Lindblad action for narrow lacunae of A, and how to predict a location of the peaks remains an intriguing 
open question. As coupling F decreases, the resonances become smoother and then ultimately disappear, but their precursors 
can still be seen for T of the order of the Heisenberg exchange interaction and larger (data not shown). 

Examining analytic expressions for the current for small N, we can extrapolate its behaviour for A 1 and < A ^ 1 and 
finite N as 

jlr^oo,A»i = ^ (a^) ^^^^ 

ilr^oo,o<A«i = 0((-l)^A^-2) for odd iV (34) 

The Eq.([33]l implies an exponential decay of the current with the system size N for fixed large A, j lp^^^ = 0(e^"^) 
with a = hi A, and is generically expected from an insulating nature of the XX Z chain in the A ::g> 1 Ising limit. Similar 
behaviour of the NESS current at large A was found for driven chain with Z-polarized boundary baths [23]. On the other hand, 
Eq.(l34li implies a complete flattening of the current near A = in the thermodynamic limit — > oo, which is a new unusual 
feature. The flattening is already well seen in Fig|4]for N = 5,7. Note that the results (|33] |. (l34l i would be very difficult to obtain 
by studying perturbative series in A or in 1/ A: in both cases A <C 1 and 1/ A <C 1 the first nonzero term will appear only in 
the order A^ — 2 of the expansion! 

C. Generic F and arbitrary driving k 

The analytic results derived in previous subsection concerned the limits of F — > oo and strong driving k = k' = 1. What can 
we say about generic F and k? Examining perturbative series (fTSl l for various A^ and A we notice that nonzero contributions to 
NESS current appear only at even orders of the expansion F~^" . If the NESS current j (F, A) is an analytic function of F, then 
we can expect the current to be an even function of F. In particular, for small F ^ 1 one expects the the current to behave as 
j(F ^ 1) = 0{r^), since the current must vanish at zero coupling F = 0. Our expectations are confirmed by exact expressions 
for the current as function of F, obtained for special case of isotropic XXX model, for which one can sum up the series ( fTSl l 
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Figure 5: Typical {(t;^) -profiles (filled circles) and (cr^)- profiles (open circles) outside the resonance (Panel (a), A = 1.8, A*" = 6) and inside 
the resonance (Panel (b), A = 0.5875, A'^ — 6). Frequency tripling is clearly observed in the Panel (b). Panel (c) shows dependence of the 
average spin projection (a^) in the middle of the chain n = (A^ + l)/2 on A, for A'' = 3, 5, 7 (from up to down). The curves for A'' = 3, 5 
are shifted vertically by 1 and 0.5 for a better view. 




Figure 6: NESS spin current j vs F for A^ = 3, 4 (thick and thin line, respectively) for isotropic Heisenberg model A = 1, given by ll35t.( l36t . 
Dashed lines indicate limiting values of the current for F — > cxd. 



and obtain 

i6r^(3 + r^) 

JliJliv=3,A=i - 27 + 52r2 + 12r4 ^ ' 

= 8r^(27 + 46r^ + ior4) 

J{ >\N=i,A=i 3(27 + 136r2 + 144r4 + 28r6) ■ ^ ' 

Indeed the above expressions contain only even T orders and j{T -C 1) = O(r^), see also Fig|6l Numerically, we observe 
j(r ^ 1) = 0{T'^) for other A values as well. Generically, we find the current amplitude to be a smooth function of F 
for a fixed value of the anisotropy A. For large couplings, it asymptotically approaches its limiting value jo{^) at F — > oo. 
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Interestingly, the magnetization current is not suppressed for large couplings T, as it happens for Z-polarized Lindblad baths 
ll23ll . ll2l1l due to a quantum Zeno effect. Indeed for our choice of the bath © and large coupling, the quantum Zeno effect 
"freezes" tr^ and spin components at the boundaries, while hoppings in the component, contributing to the magnetization 
current, may still occur. 

Finally, we investigated the current dependence on the boundary gradient k, which can be manipulated by a choice of the 
amplitudes of the Lindblad operators in (IDi. By Fourier law, we may expect the current to be proportional to k, for small gradients 
K, and introduce a conductivity tensor Xap as a proportionality coefficient between the current j = and the infinitesimal 

boundary gradient {5a^) /N in the spin component /3, j" = Xapi^jr^ 

With our perturbative method we construct iteratively an expansion (fTsl l for a NESS for an arbitrary amplitude of the boundary 
gradient k from (|4|i. We observe that the current is a rather nontrivial function of k, and in particular that it vanishes as or 
faster for small k. Consequently, limK-^o = so that the off-dia gon al elements of the conductivity tensor Xzx,Xzy 

vanish. Note that diagonal elements of the conductivity tensor are nonzero 1121 11 . Numerically, we find that the leading order of 
expansion of / k for small k depends on N and A in a nontrivial way. More details will be given elsewhere. 

Finally, note that our model shows a behaviour fundamentally different from that with usually considered boundary driving 
along the Z-axis, referred to below as a "scalar" setup. Just to mention two crucial differences: the scalar setup leads to a 
ballistic current in the critical gapless region of the XXZ model A < 1 ll23ll . while in the XY-dnven chain the alternating 
signs phenomenon ( [TOlfTSl l rules out the possibility of a ballistic current. In the "scalar" setup the current at strong coupling is 
suppressed due to a Zeno effect, j^~^oo 0, while in the XY-dn\en chain of finite size the spin current in the limit F ^ co 
remains finite. 



V. CONCLUSIONS 

We considered an open XXZ spin 1/2 chain with a twisting in XY- plane, imposed by boundary gradients. Symmetries 
of the Lindblad equation were found which impose drastic restrictions on the steady state reduced density matrix, and various 
observables. For our boundary setup the NESS energy current vanishes, and the magnetization current alternates its sign with 
system size, which rules out the possibility of a ballistic current. We argued that the current sign alternation under trasverse 
gradients is a robust phenomenon in quantum transport which does not rely on integrability. The qualitative difference in 
quantum transport in even- and odd-sized systems adds to a list of similar phenomena observed in quantum diffusion in periodic 
chains i43ll . To support our findings, we constructed a nonequilibrium steady state solution of the Lindblad Master equation in 
the form of a perturbative expansion in orders of F^^, and calculated the zeroth-order term analytically for small system sizes, 
by solving a set of linear equations which guarantee the self-consistency of the expansion, the secular conditions (ISTT l. We find 
that, for A ^ 0, the current remains finite even in the limit of infinitely strong effective coupling. 

Further on, for large couplings we find a nontrivial dependence of the magnetization current on the anysotropy i(A), and 
observed sharp peaks in j(A) inside the critical gapless region of the XXZ model A < 1, the origin of which was attributed 
to an appearance of higher Fourier harmonics in magnetization profiles. We observe an anomalous flattening of the current near 
A = point for odd N and noncommutativity of Umits A — !> and F^^ — > for even N. For weak driving k 1 and small 
system sizes, we find \imK^o j{K)/K = for all values of A, signalizing subdiffusive current in the direction trasverse to XY 
boundary gradients. 

Acknowledgements The author thanks N. Plakida, R. Vaja, M. G. Pini, M. Salerno and R. Moessner for stimulating discus- 
sions. Support from the Italian MlUR through PRIN 20083C8XFZ initiative is acknowledged. 

Appendix A: Inverse of the Lindblad dissipators and secular conditions. 

Cl and Cji are linear super-operators acting on a matrix p as defined by (|3]l. In our case, each super-operator act locally 
on a single qubit only. We find the eigen-basis {(f)%}a=i of \CR(f)'^ = \a<l)% in the form = {p/j, cr^, cr^, cr^}, where 
CT^, (T^, CT^ are Pauli matrices and p/f = ^ (/ + kct^). The respective eigenvalues are {Aq} = {0, 2, 1, 1}. Analogously, the 
eigen-basis and eigenvalues of the eigenproblem ^C^cj)^ — Hfscj)^ are 0l ~ {p^, cr^', ct^, cr^} and {^.p} = {0, 1, 2, 1}, where 
PL = 5 (-f ~ Kcr^)- Since the basises (pn and 0^ are complete, any matrix F acting in the appropriate Hilbert space can be 
expanded as 

4 4 
a=l 13=1 

in the unique way. Indeed, letus introduce complementary basises ipu = {I, ^{a^ — kI), iu*', ■^cr^} and tp^ = {I, \<^^ ^ \{'^^ ~^ 
kI), ^cr^}, trace-orthonormal to the (f>]i, 0^ respectively, Tr{tp]^(l)'^) = Sa-y, Tr{^p'J^(j)^) = Sfsj. Then, the coefficients of the 
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expansion dATJ are given by i^^Q = Tri^Ar((T/)^ (g) On the other hand, in terms of the expansion dATT ) 
the superoperator inverse [C]^ + Cji)^^ is simply 

-2{Cl + CrY^F ^Y^—^^ (l)l®Fp^®d^% (A2) 

Note however that the above sum contains a singular term with a = /3 = 1, because Ai + /^i = 0. To eliminate the singularity, 
one must require Fn = Tri^^F — 0, which produces the secular conditions ^A\ . 



[1] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, (2002). 
[2] M.B. Plenio and PL Knight, Rev. Mod. Phys. 70, 101 (1998). 

[3] F. Heidrich-Meisner, A. Honecker, and W. Brenig, Eur. Phys. J. Special Topics 151, 135 (2007), and references therein. 

[4] B. S. Shastry and B. Sutherland, Phys. Rev. Lett. 65, 243 (1990). 

[5] X. Zotos, J. Phys. Soc. Jpn. Supp. 74, 173 (2005) and references therein. 

[6] A. Klumper, Lect. Notes Phys. 645, 349 (2004). 

[7] X. Zotos, Phys. Rev. Lett. 82, 1764 (1999); J. Benz, T. Fukui, A. Klumper and C. Scheeren, J. Phys. Soc. Jpn. Suppl. 74, 181 (2005). 
[8] F. M. D. Haldane, J. Phys. C 14, 2585(1981); A. O. Gogolin and N. V. Prokof 'ev, Phys. Rev. B 50, 4921 (1994). 

[9] H. Schulz, Correlated Fermions and Transport in Mesoscopic systems, edited by T. Martin, G. Montambaux, J. Tran Thanh Van (Fron- 
tieres, Gif-sur-Yvette), 1996. 

[10] T. Antal, Z. Racz and L. Sasvari, Phys. Rev. Lett. 78, 167 (1997); T. Antal, Z. Racz, A Rakos, and G. M. Schutz, Phys. Rev. E. 57, 5184 
(1998). 

[1 1] X. Zotos and R Prelovsek, Phys. Rev. B 53, 983 (1996). 

[12] B. N. Narozhny, A. J. MiUis, and N. Andrei, Phys. Rev. B 58, R2921 (1998). 

[13] P Prelovsek, S. El Shawish, X. Zotos, and M. Long , Phys. Rev. B 70, 205129 (2004). 

[14] J.V. Alvarez and C. Gros, Phys. Rev. Lett. 88, 077203 (2002). 

[15] U Schollwock, Rev. Mod. Phys. 77, 259 (2005). 

[16] D. Gobert, C. Kollath, U. Schollwock, G.M. Schutz, Phys. Rev. E 71, 036102 (2005). 

[17] T. Prosen and M. Znidaric, J. Stat. Mech. (2009) P02035; M. Znidaric, J. of Stat. Mech. P12008 (201 1). 

[18] N. Hlubek, P Ribeiro, R. Saint-Martin, A. Revcolevschi, G. Roth, G. Behr, B. Buchner, and C. Hess, Phys. Rev. B 81, 020405 (2010). 

[19] T. Prosen, New. J. Phys. 10, 043026 (2008). 

[20] B. Zunkovic and T. Prosen, J. of Stat. Mech. P08016 (2010). 

[21] G. Benenti, G. Casati, T. Prosen, D. Rossini and M. Znidaric, Phys. Rev. B 80, 035110 (2009). 

[22] T. Prosen and M. Znidaric, Phys. Rev. Lett. 105, 060603 (2010). 

[23] T. Prosen, Phys. Rev. Lett. 107, 137201 (2011). 

[24] S. Jesenko and M. Znidaric, Phys Rev. B 84, 174438 (201 1). 

[25] X. Zotos, F Naef, and P Prelovsek, Phys Rev. B 55, 11029 (1997). 

[26] M. Znidaric, Phys. Rev.Lett. 106, 220601 (2011). 

[27] M. Znidaric, J. Stat. Mech. L05002 (2010). 

[28] M. Znidaric, J. Phys. A, 43 415004 (2010). 

[29] S. R. Clark, J. Prior, M. J. Hartmann, D. Jaksch and M. B. Plenio, New J. of Phys. 12, 025005(2010). 
[30] T. Prosen and B. Zunkovic, New J. of Phys., 12, 025016 (2010). 

[31] Strictly speaking, the boundary spins become fully antiparallel in the limit of infinitely large coupling with reservoir and strong boundary 
driving. 

[32] V. Popkov, M. Salerno and G. M. Schutz, Phys. Rev. E 85, 031 137 (2012). 
[33] T. Prosen, Phys. Rev. Lett. 106, 217206 (2011). 

[34] D. Karevski and T. Platini, Phys. Rev. Lett. 102, 207207 (2009); S. Attal and Y. Pautrat, Ann. Henry Poincare 7, 59 (2006). 

[35] Within the repeated interaction scheme (29ll . (34ll . the LME l[T) describes a system evolution where coupling to reservoir is rescaled as 

r / -^/r where r is a time interval between the two successive interactions with the ancilla bath. 
[36] V. Gorini, A. Kossakowski and E. C. G. Sudarshan, J. of Math. Phys. 17, 821 (1976). 

[37] The relaxation times depend on spin component, the targeted spin component relaxing two times faster than other spin components. E.g. 

at the left boundary, x- and 2:-spin components of p relax as e^^^ while the y-sp'm component relax as e^^'^*, see 13211 . 
[38] If K = ft' < 1, the situation is more subtle since then the target states pL and pR are mixed 
[39] D.E. Evans, Comm. Math. Phys, 54, 293 (1977) 
[40] T. Prosen, Phys. Scr. Lett. 86, 058511 (2012). 
[41] A. Heinrich and S. Loth, Science 332, 1039 (2011). 

[42] M. Suzuki, Physica 51, 277 (1971) 28; V. Aksenov, M. Bobeth, N. Plakida, J. Schreiber, J. Phys. C. 20, 375 (1987). 
[43] A. C. de la Torre, H. O. Martin and D. Goyeneche Phys. Rev. E 68, 03 1 103 (2003). 



